#include<iostream>
#include<string>
#include<cmath>
#include<vector>
#include<cstdlib>
#include<algorithm>
#include"mpi.h"

using namespace std;

const int N = 2000;
int mid, proc;
void copyMemory(const double *M, double *A) 
{
	for (int i = 0; i < N*(N + 1); i++)
		A[((i / (N + 1)) % proc)*( (N * (N + 1)) / proc) + (N+1) * ((i / (N + 1)) / proc) + i % (N + 1)] = M[i];
}
int main(int argc,char* argv[]) {
	MPI_Init(&argc,&argv);
	double t1, t2;
	t1 = MPI_Wtime();
	MPI_Comm_rank(MPI_COMM_WORLD,&mid);
	MPI_Comm_size(MPI_COMM_WORLD,&proc);
	double *M, *c, *M_trans,*A,*tempRow;
	int *map = new int[N]; int *rmap = new int[N];
	for (int i = 0; i < N; i++) {
		map[i] = -1; rmap[i] = -1;
	}
	if (mid == 0) 
    {
		M = new double[N*(N+1)];
		M_trans = new double[N*(N + 1)];
		for (int i = 0; i < N*(N + 1); i++) {
			M[i] = rand() % 10;
		}
		
		copyMemory(M, M_trans);
		A = new double[N*(N + 1) / proc + 1];
		c = new double[N];
	}
	else {
		M = new double[1];
		M_trans = new double[1];
		A = new double[N*(N + 1) /proc+1];
		c = new double[N];
	}
	tempRow = new double[N + 1];
	MPI_Scatter(M_trans, N*(N+1) / proc,MPI_DOUBLE,A,N*(N+1)/proc,MPI_DOUBLE,0,MPI_COMM_WORLD);
	delete M_trans;
	double local_max_value, global_max_value;
	int local_max_id, max_proc_id, max_pro_id_1, global_max_id;
	for (int row = 0; row < N-1; row++) {
		local_max_value = -1e10; local_max_id = -1;  max_pro_id_1 = -1, global_max_id=-1;
		for (int i = 0; i < N / proc; i++) {
			if (rmap[i * proc + mid] < 0) {
				if (A[i*(N + 1) + row] > local_max_value) {
					local_max_value = A[i*(N + 1) + row];
					local_max_id = i;
				}
			}
		}
		MPI_Allreduce(&local_max_value, &global_max_value, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
		if (global_max_value > local_max_value-0.00001&&global_max_value < local_max_value+0.00001) 
		{
			max_pro_id_1 = mid;
			for (int i = 0; i < N + 1; i++) {
				tempRow[i] = A[local_max_id*(N + 1) + i];
			}
			global_max_id= local_max_id * proc + mid;			
		}
		MPI_Allreduce(&max_pro_id_1,&max_proc_id,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
		MPI_Bcast(&max_proc_id,1,MPI_INT,0,MPI_COMM_WORLD);
		MPI_Bcast(tempRow,N+1,MPI_DOUBLE,max_proc_id,MPI_COMM_WORLD);
		MPI_Bcast(&global_max_id,1,MPI_INT,max_proc_id,MPI_COMM_WORLD);
		
		map[row] = global_max_id; rmap[global_max_id] = row;
		for (int i = 0; i < N/proc; i++) 
		{
			if (rmap[i * proc + mid] < 0) {
				double temp = A[i*(N + 1)+row] / tempRow[row];
				for (int j = row; j < N + 1; j++) {
					A[i*(N + 1) + j] -= tempRow[j] * temp;
				}
			}
		}
	}
	
	MPI_Gather(A,N*(N+1)/proc,MPI_DOUBLE,M,N*(N+1)/proc,MPI_DOUBLE,0,MPI_COMM_WORLD);
	if (mid == 0) {
		for (int i = 0; i < N; i++) {
			if (rmap[i] == -1) {
				map[N - 1] = i;
			}
		}
		printf("\n");
	}
	if (mid == 0) 
	{
		for (int i = N-1; i >=0; i--) {
			int index = map[i] % proc*(N / proc) + map[i] / proc;
			for (int j = N-1; j >i; j--) {
				M[index*(N+1)+N] -= c[j] * M[index*(N + 1) + N - j];
			}
			c[i] = M[index*(N + 1) + N] / M[index*(N+1)+i];
		}
		for (int i = 0; i < N; i++) {
			printf("C[%d]=%.2f\n",i, c[i]);
		}
	}
	t2 = MPI_Wtime();
	printf("it take %.2lfs\n", t2-t1);
	MPI_Finalize();
}